Density profiles of a colloidal liquid at a wall under shear flow 
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Using a dynamical density functional theory we analyze the density profile of a colloidal liquid near 
a wall under shear flow. Due to the symmetries of the system considered, the naive application of 
dynamical density functional theory does not lead to a shear induced modification of the equilibrium 
density profile, which would be expected on physical grounds. By introducing a physically motivated 
dynamic mean field correction we incorporate the missing shear induced interparticle forces into the 
theory. We find that the shear fiow tends to enhance the oscillations in the density profile of hard- 
spheres at a hard-wall and, at sufficiently high shear rates, induces a nonequilibrium transition to a 
steady state characterized by planes of particles parallel to the wall. Under gravity, we find that the 
center-of-mass of the density distribution increases with shear rate, i.e., shear increases the potential 
energy of the particles. 
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I. INTRODUCTION 

Classical Density Functional Theory (DFT) provides a 
powerful and general framework for determining the equi- 
librium microstructure and thermodynamics of classical 
many particle systems [J 0] • Of particular interest is the 
one-body density profile p(r) resulting from application 
of a time- independent external potential field ^"^'(r). 
Within DFT, the density profile of a one-component sys- 
tem follows from functional minimization of the Grand 
potential 

n[p]^T,^[p]+T,^[p]+ y'dr(F-'(r)-A^)p(r), (1) 

with respect to p(j), where fi is the chemical potential 
and -Fox[p] is the unknown 'excess' part of the Helmholtz 
free energy, containing details of the interparticle inter- 
actions. The ideal gas free energy is given exactly by 



drp(r)[ln(A'V(r)) - 1], 



(2) 



where A is the thermal de Broglie wavelength. For many 
important model systems (e.g. hard-spheres [1, colloid- 
polymer mixtures rod-sphere mixtures ^) there 
exist accurate approximations for the excess Helmholtz 
free energy which yield equilibrium density profiles in 
excellent agreement with those obtained from numerical 
simulation. 

Given that DFT provides a clear method for deter- 
mining equilibrium density distributions, it is natural to 
next investigate the dynamics of the density profile in 
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the presence of a time-dependent external field y™*(r, t). 
The simplest, phenomenological, route to an equation of 
motion for p{v,t) is to assume that the average particle 
current j(r,<) arises from the gradient of a nonequilib- 
rium chemical potential 



j{r,t) = -Tp{r,t)Vpir,t), 



(3) 



where F is the mobility. Assuming that p{Y,t) is given 
by the functional derivative of the Helmholtz free en- 
ergy with respect to p(r, t) and employing the continuity 
equation thus leads to the familiar equation of dynamical 
density functional theory (DDFT) 



dp{-r,t) 
dt 



Fp(r,t)V 



6p{v,t) 



(4) 



where T is the equilibrium Helmholtz free energy func- 
tional, evaluated using the instantaneous nonequilibrium 
density profile. Although equation (|4]) was first pro- 
posed by Evans [l| , subsequent work has clarified greatly 
the nature of the approximations involved. Both Mar- 
coni and Tarazona 0, ll], and Archer and Evans [§], 
have demonstrated that approximating the nonequilib- 
rium chemical potential using the equilibrium Helmholtz 
free energy is equivalent to assuming that the inhomoge- 
neous pair correlations in nonequilibrium are the same as 
those for an equilibrium fluid with density profile p(r,t). 
Specifically, for a system interacting via a pair-potential 
0(|ri — r2|) the equilibrium sum-rule [l| 



J dr2P^^Hrur2)Vi(j){\ri-r2\) = p(ri)Vi 



6p{ri 



(5) 



is assumed to hold. The integral on the l.h.s. of ([5]) oc- 
curs when coarse graining the iV-particle Smoluchowski 
equation to the one-body level by integration over A^ — 1 
degrees of freedom. Applying the equilibrium equality 
^ to close the resulting nonequilibrium expression leads 
directly to The implicit 'adiabatic' approximation in 
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applying ([5]) to noncquilibrium is that the one-time pair 
correlations arc instantaneously equilibrated to those of 
an equilibrium system with density p{r,t). For a wide 
range of systems, the good qualitative agreement be- 
tween the results of Brownian dynamics simulation and 
DDFT validates the adiabatic approximation when ap- 
plied to inhomogeneous fluid states out of equilibrium. 
However, the approximation may break down for dense 
fluids close to dynamical arrest (e.g. hard-sphere-like col- 
loidal glasses), for which the structural relaxation time 
becomes large. 

The possibility of going beyond the adiabatic approxi- 
mation has been explored on a formal level [l^. However, 
an explicit and implementable method of incorporating 
temporal nonlocality into the theory remains to be found. 
More recently, the DDFT (jU has been rederived using 
projection operator method s 11 ll a nd extended to incor- 
porate pair hydrodynamics jl3l - ll5| . orientational degrees 
of freedom [Tfl] and even self-propelled particles • Go- 
ing beyond overdamped Smoluchowski dynamics, Mar- 
coni and coworkers have developed a DDFT-type ap- 
proach to treat inertial systems which makes possible the 
study of thermophorcsis [12 1. 

In order to address the influence of external flow on 
inhomogeneous density distributions, a DDFT incorpo- 
rating the advection of particles by a flowing solvent 
has been developed [l^. The resulting advcctcd-DDFT 
equation of motion has a form identical to that of the 
standard theory ([4]), but with the time derivative re- 
placed by the Stokes derivative 



dt 



+ V-(p(r,t)v(r,i)) 



rp(r,i)V 



6p{v,t) 



(6) 



where v(r, t) is the velocity of the solvent. The advected- 
DDFT is therefore simply standard DDFT in the comov- 
ing frame and is thus subject to the same adiabatic ap- 
proximation as the original theory. However, compared 
to situations for which the relaxation dynamics is of a 
purely diffusive nature, application of the equilibrium 
identity ([S]) to an externally driven system represents 
a more severe approximation. For example, in the ab- 
sence of an external potential field, equation ([5|) is auto- 
matically, and trivially, satisfied for a homogeneous and 
isotropic fluid. This is not the case for a driven system. 
Even when y°'^'(r) = 0, the presence of an external flow 
field distorts the pair correlation functions and renders 
the integral term finite, whereas the spatial homogene- 
ity of the one-body density yields a value of zero for the 
r.h.s. 

There is considerable current research activity in the 
development of theoretical approaches to treat colloidal 
systems driven into nonequilibrium states by external 
flow. Much of the focus has been on the description 
of dense states close to the glass transition (see e.g. 
[l^H^]). While the mode-coupling- type approaches em- 
ployed in these studies are capable of capturing noner- 



godic behaviour, their application is restricted to sys- 
tems with a spatially homogeneous density distribution. 
In contrast, Eq.® is ideally suited to address spatial 
inhomogeneity, but is incapable of describing glassy dy- 
namics. 

In the present work we will consider the application of 
^ to driven steady states. In order to clearly expose 
the nature of the underlying approximations we will fo- 
cus on the specific test case of interacting colloidal parti- 
cles at a planar wall, with a shear flow acting parallel to 
the wall. Consideration of this particular external field 
and flow geometry reveals a serious deficiency of apply- 
ing (O to close the equation of motion for the one-body 
density of driven states. The physics which is lost in 
making the closure approximation arises from a coupling 
between the intcrparticlc interactions and the external 
flow field and would, in an exact treatment, be captured 
implicitly by the flow induced distortion of p'^^)(ri, r2). 
Previous studies based on © have focused on two cases: 
(i) Noninteracting particles, for which the only relevant 
coupling is that between the external potential and flow 
fields [2l|, [23 , (ii) Spherically inhomogeneous soft Gaus- 
sian particles [13, Hj, [13 • In these studies the combina- 
tion of soft particle interactions and spherical inhomo- 
geneity served to obscure the failings of Eq.® addressed 
in the present work. 

The paper will be structured as follows: In section [H] 
we specify the system under consideration. In section Hill 
we introduce the problem presented by the absence of a 
coupling between the external flow and the interparticle 
interactions. In sections IIVI and |V] we develop a mean 
field theory which captures the desired coupling and pro- 
pose a simple approximation for the required convolution 
kernel. In section IVll we detail the Rosenfeld functional 
used to approximate the excess free energy functional. In 
section I VIII we present and analyze the density profiles 
of hard-spheres at a hard wall, both in the presence and 
absence of gravity. In section IVIIII we give a discussion 
and provide an outlook for future work. 

II. THE SYSTEM 

We consider a system of N spherical colloidal parti- 
cles dispersed in an incompressible Newtonian solvent at 
temperature T . The diameter d of the strongly repulsive 
colloidal core provides the basic unit of length. Assuming 
that the colloidal momenta arc instantaneously thermal- 
ized, the time evolution of the probability distribution of 
particle positions, ^'(t) = ^'({ri},t), is dictated by the 
Smoluchowski equation [25j 



dMt)_ 
dt 



= 0, 



where the probability flux of particle i is given by 



(7) 



(8) 



3 



o9 o 



a- 
o 



(a) 



>Soooo o" 



neo o OOP g 



(bL 



^ 


o 


(c) 



FIG. 1: A schematic two-dimensional illustration showing 
the effect of shear flow on the microstructure of hard spheres 
at a hard wall, (a) A typical equilibrium configuration, (b) 
Shear flow leads to the formation of layers in the a^z-plane 
as particles at different values of y attempt to overtake each 
other, (c) Focusing on a binary collision close to the wall. 
The particle closer to the wall gets pushed against it, while 
the colliding particle is forced to 'roll around' the other in 
order to move past. 



with l3 = l/fc^T. The hydrodynamic velocity of particle 
i due to the applied flow is denoted by Vi(i) and the 
diffusion tensor Dij describes (via the mobility tensor 
Tij ~ PDij) the hydrodynamic mobility of particle i 
resulting from a force on particle j. The force Fj on 
particle j is generated from the total potential energy 
according to = —djUN and includes the influence of 
an external one-body potential field, as well as the forces 
generated by interaction with other particles (taken here 
to be pairwise additive) 

C/^({rJ,t)^^l^-*(r„0 + E'^(|r^'-r^l)- (9) 

i i<j 

The three terms contributing to the flux thus represent 
the competing effects of (from left to right in ([5])) external 
flow, diffusion and potential interactions. 

In order to focus on the thermodynamic (as opposed 
to hydrodynamic) aspects of the cooperative particle mo- 
tion we will neglect hydrodynamic interactions (HI) be- 
tween the colloids. The expression ([5]) is thus approxi- 
mated in two ways: (i) The influence of the TV-particle 
configuration on the mobility of a given particle is ne- 
glected, Dij = DoSij, where Dq is the bare diffusivity. 
(ii) The velocity field is determined by the translationally 
invariant (traceless) velocity gradient tensor describing 
the affine motion, = v{ri,t) = /t(t) • r^. Neglecting 
the dependence of K{t) upon the colloid configuration 
enables the externally applied flow to be prescribed from 
the outset, without requiring that this be determined as 
part of a self consistent calculation. 



III. INTERACTION-FLOW COUPLING 

In the present work we wish to study the influence 
of flow on a fully interacting, inhomogeneous system. 
So far, the only application of (O has been to spher- 
ically inhoi TLOg eneous systems subject to a variety of 
flow fields |l8l . [2l| - |23 | (representing a single colloid in 
a sea of ideal or Gaussian polymers). In particular, 
under shear flow, the ideal particles accumulate at the 



colloid surface within the two compressive quadrants 
{I = r-{K + k'^) • r < 0) and are depleted from the exten- 
sional quadrants (7 > 0), leading to 'wake' formation at 
larger shear rates [22] . The advected DDFT ^ thus cap- 
tures a certain coupling between inhomogeneities in the 
density field and external flow. However, this 'external 
potential-flow coupling' is a straightforward consequence 
of employing an external potential (e.g. a fixed sphere) 
which either directly obstructs the particles as they at- 
tempt to follow the affine flow, or perturbs the solvent 
flow field such that the particles are swept around the 
obstacle (when the appropriate Stokes flow is employed). 

A more demanding and illuminating test of the 
advected-DDFT approach is provided by considering ex- 
ternal potentials which do not directly hinder the affine 
path of the advected particles, but which may neverthe- 
less be expected to generate flow dependent density pro- 
files. Emphasis may thus be placed upon the nontriv- 
ial 'interaction-fiow coupling'. For the purpose of the 
present work we will thus focus on the special case of 
a time-independent external potential which varies in y- 
direction and restricts the particles to move in the half 
space 



y < d/2 
y > d/2. 



(10) 



The translational invariance of V^^*' within the a;z-plane 
has the consequence that the resulting equilibrium den- 
sity distribution varies in y-direction only. In addition 
to the external potential field, we specialize the external 
flow field to a steady simple shear with flow in a;-direction 
and a linear gradient in the y-direction 



yje^ 



(11) 



with shear rate 7 (corresponding to a velocity gradient 
tensor Kap = jSaxSiSy). The pair potential entering ([9]) 
represents a hard-sphere repulsion 




r < d 
r > d. 



(12) 



The situation under consideration is shown schematically 
in the second panel of Fig[TJ Note that the zero-shear 
plane may be located at ?/ = without loss of generality, 
as only relative particle velocities are physically relevant. 

Application of ^ to treat the specified nonequilibrium 
situation immediately reveals the problem at hand. We 
have already noted that our chosen geometry leads to 
translational invariance of the equilibrium density distri- 
bution within the a;z-plane, poq(r,t) = Pcq{y,t). For the 
shear flow (fTT|) the advective term in ^ is thus given by 

V • (p(r, t)v(r, i)) = V • (j/ 7 Piy. i) e.) = 0. (13) 

Within the advected-DDFT approach the applied shear 
flow has no influence on the density profile at the wall. 
Equation ^ thus reduces to which, for the present 
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time-indcpcndcnt external potential, yields the equilib- 
rium density profile. This disappointing conclusion con- 
tradicts physical intuition and presents a fundamental 
failing of Eq.®. 

In the right panel of FigH] we sketch what we con- 
sider to be the correct physical picture. As a shear flow 
is applied parallel to the wall the particles experience a 
(nonconservative) force proportional to their perpendicu- 
lar separation from the wall. Particles at different values 
of y thus seek to move past each other, perturbing the 
equilibrium microstructure and leading, at higher shear 
rates, to the formation of particle layers in the xz-plane. 
In particular, when a pair of particles collide close to the 
wall, the particle at smaller y will be forced against the 
wall, whereas the second particle will be forced to roll 
around its neighbour in order to follow as closely as pos- 
sible the affine solvent flow. These physical mechanisms 
should be manifest in the nonequilibrium density profile. 

Wc note that Brownian dynamics simulations |39l | dis- 
play two dimensional particle layering at intermediate 
shear rates, followed by an additional symmetry break- 
ing in the vorticity direction at high shear rates, charac- 
terized by the formation of particle chains in x-direction. 
It is important to realize that the density distribution 
described by DFT represents an average over all z- 
coordinates of the particle chains (which, for an infinite 
system, are not pinned to a particular location in z). Our 
assumption that p{r, t) = p{y, t) is thus fully justified for 
the present density functional based study. 



IV. MEAN FIELD THEORY 

In a system without HI, A^-particle configurations for 
which colloidal particles overlap are forbidden by the in- 
finitely repulsive colloidal pair potential. While an exact 
statistical mechanical treatment (i.e. exact evaluation of 
the integral in ([5])) would lend such unphysical configu- 
rations zero statistical weight, care must be exercised in 
approximate treatments which do not fully satisfy this 
important geometrical constraint. In the present situa- 
tion it would appear that application of the equilibrium 
sum rule ([5]) does not satisfy exactly the no-overlap 'core 
condition' when applied to driven states. However, it is 
not at all clear how to improve the approximation ([5]) 
in a way that both incorporates the flow induced distor- 
tion of p^^^(ri, r2, t) and enables its weighted integral to 
be approximated using an equilibrium free energy func- 
tional. For this reason we take an alternative route and 
attempt to include the missing physics by modifying the 
advective term in This approach is motivated by 

considering the dynamics of hydrodynamically interact- 
ing dispersions. 

In a system with HI, the hydrodynamic velocity en- 
tering ([8]) can be decomposed into affine and particle in- 
duced fiuctuation terms Vi(<) = nit) ■ -f- vf (t), where 
vf (t) can be expressed in terms of the third rank hy- 
drodynamic resistance tensor [26j . The fluctuation term 



describes the disturbance of the affine solvent flow by the 
particles and ensures that a pair of approaching particles 
'flow around' each other, without coming into direct con- 
tact. The impenetrable character of the particles, repre- 
sented by the pair potential (|12p . thus enters indirectly 
by providing a boundary condition for the solvent flow. 
Integration of ([7]) over A^— 1 degrees of freedom yields an 
advective term 



V ■ p{v,t)[K{t) 



v«(r,i): 



(14) 



where v*'(r,t) = (vf(t)) is the conditional average over 
N — \ coordinates, given that a particle is located at r. 
In contrast to p^. the divergence in dU]) is not necces- 
sarily zero for the inhomgeneous system presently under 
consideration. For hydrodynamically interacting systems 
the fluctuation term may thus provide the desired con- 
tribution to the flux in y-direction. 

The above considerations suggest a possible solution to 
the problem posed by (|13p . By empirically incorporating 
a conditionally averaged fluctuation term into the veloc- 
ity fleld driving our Brownian system, we aim to mimic 
the hydrodynamic fluctuation term discussed above. In 
this way we can correct, at least to some extent, the oc- 
curance of unphysical overlaps by enforcing that there be 
no radial flux between pairs of particles at contact. We 
envisage a system of hard-core particles under shear flow 
in which there occur frequent and random binary colli- 
sions. At each collision the particles rotate around each 
other according to some specified rule (for our specific 
choice, see section |V] below) before moving apart along 
their respective streamlines. For a homogeneous system, 
this mechanism gives rise to zero net flux through any 
given plane at constant y. However, in the presence of 
the external potential field pU)) , the inhomogeneous den- 
sity distribution will lead to a y-dependent flux which will 
modify the equilibrium distribution, until it is balanced 
by an equal and opposite diffusive flux, thus establishing 
a nonequilibrium steady state. 

As the density proflle under consideration varies only 
in y-direction, we need only seek the y-component of the 
fluctuation contribution. The arguments presented above 
thus suggest the mean field term 



vl{y.t) 



dy'p{y\t)t{y-y') 



(15) 



which expresses the contribution of the density at y' to 
the average velocity at y. The 'flow kernel' Vy{y) entering 
([T5|) describes the y-component of the velocity of a parti- 
cle which collides with a neighbour at vertical separation 
y. Our modified version of (jS]) thus becomes 



dp{y,t) 
dt 



d_ 

dy 



p{y,t)vl{y,t) 



+ Pe-^p{y,t) 



d ST[p{y,t)] 
dy Sp{y,t) 



,(16) 
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where we have scaled length using d and time using 7, 
leading to an explicit appearance of the Peclet number, 
Pe = jcP / Dq, expressing the competition between affine 
advection and diffusive motion. The steady state solution 
of (fT6)) is given by 



p{y) = zexp[-/3V''''\y)+c 



(1) 



{y) + Pe dy'vliy'] 
•'v 



(17) 



where z is the equilibrium fugacity and we have intro- 
duced the one-body direct correlation function [l[ 



.(1) 



(r) 



5p{y,t) 



(18) 



The integral in (|17|) can be regarded as a noncquilibrium 
contribution to the intrinsic chemical potential. It should 
be noted that the requirement that a homogeneous den- 
sity generates no particle flux in y-direction implies that 
the spatial integral over Vy{y) is zero. For homogeneous 
states Vy{y, t) is thus zero and the bulk chemical potential 
is unchanged from that in equilibrium. 



V. FLOW KERNEL 

In order to derive the flow kernel in Eq. ()15p . we con- 
sider the relative velocities of two particles, a tagged par- 
ticle with position Yt and a reference particle at r,., dur- 
ing a scattering event, with Ar = — r,.. As discussed 
above, we neglect hydrodynamic interactions to keep the 
description as simple as possible. The incorporation of 
hydrodynamic interactions should in principle be possi- 
ble in our approach. The following derivation of the scat- 
tering velocities is based on the assumption that in any 
situation, the particles adjust their velocities such that 
they minimize the friction with the solvent. The parti- 
cles are at contact while they move around each other 
and then pass each other (see the right panel of FigH]). 
During the contact period, they have a nonvanishing ve- 
locity relative to the solvent. We will now derive the 
velocity v'^ of the tagged particle in the frame comoving 
with the reference particle, which is assumed to move 
with constant velocity 7j/r, i-e., we keep yr fixed. In a 
real scattering event one particle would move up and the 
other down. In our mean field picture, where the tagged 
particle moves in the fixed density distribution of other 
particles, we have to keep the y-position of the reference 
particle fixed. The velocity with minimal friction follows 
from minimization of the velocity relative to the solvent 
velocity v, which we denote Av, 



(v^■-v)2 = (^,^^-7AJ/)2 + ^;f +„; 



,fe2 



(19) 



As long as the particles are in contact, they move on a 
trajectory with constant distance from each other, lead- 
ing to 



Ar^ 



(20) 



Differentiation of this expression and elimination of Ace 
leads to 



, w/;Ay -I- w^Az 

V = 



^d^ ~ A?/2 - Az2 

Inserting in Eq. (fT9|) yields 



(21) 



(22) 

We can now minimize this expression with respect to 



and w^- This yields. 



My 



i;^(A2/,Az) = -7(^^j ^ d^ - ^y^ - (23) 

In order to average this value over all possible values of 
Az, we substitute 



Az = \J d^ — Aj/2 sin < 



(24) 



into ([23| and integrate over to obtain (the factor of ^ 
is inserted as normalization, and the range of the integral 
is restricted to the half circle since the velocity is zero for 
101 > Tr/2) 



(25) 



This result does not account for the density of particles 
at contact, relative to that in bulk. In order to include 
information about the local microstructure about a ref- 
erence particle, the flow distorted inhomogeneous pair 
distribution function ^(r, y) should, in principle, be in- 
cluded as a prefactor in before angular integration. 
Given that g{r,y) is not known, we make the zeroth or- 
der approximation g{v,y) = geq{r), to arrive at our final 
result 



v'yiAy) 
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--9cq{d) 



f 1 ^ 



Ay2 



(26) 



For hard-spheres the Carnahan-Starling expression 
9eq{d) = (l — 0/2)/(l — 0)'^ provides a simple and accurate 
expression for the contact value . For other choices of 
interaction potential gcqid) may be obtained using either 
integral equation methods [28| or, more consistently, an 
equilibrium test-particle calculation employing the same 
Helmholtz free energy as that used to generate the dy- 
namics. 



VI. EXCESS FREE ENERGY 



Given the equation of motion ([T6| and fiow kernel <\26^ , 
we need to specify a particular approximation to the ex- 
cess free energy functional in order to arrive at a closed 



theory for the density profile. For hard-sphere fluids 
the Rosenfeld functional 3 yields accurate results for 
both the microstructure and thermodynamics. Within 
the Rosenfeld approximation the excess free energy of 
hard-spheres is given by 



•^cxW = -noln(l-n3) + 



nin2 -ni ■n2 



1 - ns 



+ 



«2 - 3n2(n2 ■ n2) 
247r(l-n3)2 ' 



(27) 



where the weighted densities are given by convolutions of 
the density profile 



(28) 



The weight functions are characteristic of the geometry 
of the particles 





= e(r-i?), 




= Sir-R), 


c^(i)(r) 


S{r - R) 


AttR ' 


u;(o)(r) 


5{r - R) 


4ttR^ ' 


a;(2)(r) 


= -S(r~R), 
r 


a;(i)(r) 


rS{r - R) 


r inR 



(29) 



where R = d/2 is the sphere radius. Although improved 
versions of the Rosenfeld functional do exist [231, the 
original version 0] will prove sufficient for the present 
application. 



VII. RESULTS 



A. Hard-spheres 
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FIG. 2: Steady state density profiles of a hard-sphere fluid 
at volume fraction (j> = 0.45 for Pe = 0,3.142,6.283,9.425 
and 12.566. As the Peclet number is increased the oscillatory 
structure of the profile becomes more pronounced, reflecting 
the formation of particle layers in the a;2-plane. Inset (a) 
focuses on the region close to the wall, where the contact value 
p{d/2) increases linearly with Pe. Inset (b) shows the density 
dependence of the coefficient determining the nonequilibrium 
contribution to the reduced osmotic pressure /?nnc(0,-Pe) = 
'•^Pe, as determined from the contact value. 



value (see inset (a)) and height of the oscillatory peaks, 
which is accompanied by an increase in the depth of the 
minima. The enhanced structure of the profile is a direct 
consequence of the collision mechanism built into the ad- 
vective term of our theory (c.f. FigHJc) and indicates the 
development of particle layers in nonequilibrium steady 
states at finite Pe. Despite the highly structured charac- 
ter of the nonequilibrium profiles, it should be noted that 
the adsorption (i.e. the spatial integral of p{y) — pb) re- 
mains independent of Pe, where pb = is the bulk col- 
loid density. While this is a straightforward consequence 
of the continuity equation underlying ^ , it nevertheless 
provides a useful check for our numerical results. 



We first address the problem of hard-spheres at a hard 
wall {V^^^{y) = 0). The steady state equation ([TT]) was 
solved numerically (Picard iteration) using the flow ker- 
nel ((26)) and the Rosenfeld approximation to the excess 
free energy. The contact value of the radial distribution 
function employed in ()26|) was taken from the Carnahan- 
Starling equation of state [l^l . 

1. Intermediate shear rates 

In FigI2]we show density profiles calculated for a vol- 
ume fraction (p = 0-45 at various (low to intermediate) 
values of the Peclet number. In equilibrium, Pe = 0, 
the density profile shows a typical oscillatory structure 
arising from local packing contraints at the wall. Apply- 
ing a finite Pe leads to an increase in both the contact 



2. Osmotic pressure 

For hard-spheres at a hard-wall, the contact value of 
the density profile satisfies the sum rule 

/3n = p(d/2), (30) 

where 11 is the osmotic pressure. While the sum rule 
pop is generally applied to equilibrium, there seems to 
be no reason why it should not be equally valid for the 
present nonequilibrium situation (although, as far as we 
are aware, there currently exists no mathematical proof 
of this assertion). Our numerical calculations show that, 
for a given volume fraction, the contact value p(d/2) in- 
creases linearly over a the entire range of Peclet numbers 
investigated {Pe = — > 10 for each volume fraction con- 
sidered). Employing the sum rule pOp we thus find that 
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FIG. 3: Steady state density profiles for volume fraction 
(j> = 0.42 at Pe = 16.965, 17.279, 17.593, 17.907 and 18.221. 
For clarity the profiles have been translated vertically. As Pe 
is increased the profiles display an increasingly slow oscillatory 
decay to the bulk density. At a critical value of Pe = Pccrit 
the Brownian motion is no longer able to restore the equilib- 
rium structure far from the wall and shear effects dominate. 
For Pe > Pecrit the oscillations no longer decay and the entire 
sample enters a layered state, characterized by a well defined 
oscillation amplitude. The width of the laning oscillations 
away from the wall (and which extend throughout the en- 
tire sample) defines an order parameter W characterizing the 
nonequilibrium transition. The inset focuses on a single peak 
{(f) = 0.42) within the layering region for Pe = 17.907 and 
18.221. The layering peaks can be well approximated by a 
Gaussian. 

the numerically obtained osmotic pressure obeys the fol- 
lowing relation 

/3n((/), Pe) - /3neq((/)) + a{4')4>''Pe (31) 

where Q!(0) is a volume fraction dependent coefficient. 

The definition of q;(0) in the second term of (|3ip is mo- 
tivated by the exact low volume fraction results of Brady 
and Morris [s^. By solving exactly the pair Smolu- 
chowski equation in the low density limit for hard-spheres 
without HI it has been shown that the osmotic pressure 
(obtained from the trace of the stress tensor) is given by 

/?n((/) ^ 0, Pe) = pfc + A'^'^e (32) 

In inset (b) of Figl5] we show the volume fraction de- 
pendence of a. Gratifyingly, the fact that a exhibits a 
low volume fraction plateau confirms that the present 
theory indeed captures the correct scaling (~ cfP'Pe) of 
the flow induced correction to the osmotic pressure. The 



fact that we recover the correct low density scaling is a 
nontrivial output of our approach. However, the limiting 
value a{4) — > 0) = 0.164 obtained in the present work dif- 
fers from the exact value of 4/37r^ = 0.135 by a factor of 
1.2. Given the rather severe approximations employed in 
the present work, namely the mean-field term (|15p and 
flow kernel (|26l). it should not be surprising that there 
is some deviation from the exact result. Nevertheless, 
the recovery of the correct low density scaling is reas- 
suring and suggests that performing calculations with a 
renormalized Peclet number Pe* = Pe/1.2 may be ap- 
propriate, should a detailed comparison with simulation 
or experiment be required. 

We note that time-dependent solutions of (|16p (not 
considered in the present work) would enable study of 
the transient behaviour of the Osmotic pressure resulting 
from time-dependent changes in the applied shear flow, 
e.g. the onset of steady shear flow (sst . 



3. Laning transition 

Turning now to higher values of Pe, we show in Fig [3] 
density proflles for (j) = 0.42 and Pe = 16.695 18.221, 
focusing on the layering structure away from the direct 
vicinity of the wall (the region j/ > 5 is shown). As Pe is 
increased from zero to values around 17.6 the oscillatory 
structure shows an increasingly slow decay with distance 
from the wall, indicating that Brownian motion is gradu- 
ally succumbing to the influence of the applied shear flow. 
For Pe > 17.6 the decay length of the oscillatory pro- 
flles shows a strong sensitivity to variations in the Peclet 
number and diverges at a critical value Pccrit ~ 17.8. 
This divergence signifles the onset of an ordered phase 
for which the asymptotic density proflle is characterized 
by a well deflned period and amplitude of oscillation. 

The emergence of an infinitely extended oscillatory 
profile at a critical value of the Peclet number is a non- 
trivial prediction of the present theory and signifies a 
non equilibrium transition to an inhomogeneous steady 
state. Such layered states have been observed in Brow- 
nian dynamics simulations [36j but have thus far re- 
mained inaccessable to microscopically based theories. 
For Pe > Pecrit it is of interest to look at the structure 
of the individual oscillations within the layered phase. 
In the inset to Figl3] we show a single density peak at 
y w 51.4 for Pe = 17.907 and 18.221. For larger values 
of Pe the peak becomes both narrower and higher, re- 
flecting the reduced influence of Brownian motion, which 
acts to damp the oscillations and restore equilibrium. 
The density peaks in the layering region may be well 
approximated by a Gaussian, implying the existence of 
two-dimensional particle planes at high Pe values, with 
harmonic restoring forces acting against random out-of- 
plane fluctuations. 

Shear induced layering phases, similar to those pre- 
dicted by the present theory, have been observed in both 
colloidal experiments |3ll435l | and Brownian dynamics 



simulations of hard-sphere dispersions [36[. More re- 
cently, experiments on noncolloidal dispersions (no Brow- 
nian motion) under oscillatory shear have shown that the 
presence of irreversible processes when the particles col- 
lide can give rise to self-organization and the formation 
of particle lanes [stI ]. 



4- Phase diagram 

The oscillation amplitude of the density in the limit 
y ^ oo serves as an order parameter characterizing the 
nonequilibrium transition from a locally layered state, 
homogeneous in bulk, to a fully macroscopic layered 
phase. Specifically, W = max(p(?/ oo)) — min(p(?/ ^• 
oo)) provides a suitable order parameter (see Fig|3]). In 
Fig. 4 we show the nonequilibrium phase diagram in the 
{(j) , Pe) plane, obtained by examination of as a func- 
tion of Pe. For each volume fraction density profiles were 
calculated on a large grid extending beyond 300 particle 
diameters. For Pe < Pccrit the converged profiles clearly 
decay to zero as a function of y, well within our sam- 
ple size (as for the profiles for Pe = 16.965, 17.279 and 
17.593 in FigO. For Pe > Pe„it iteration of Eq.p7|) 
results in a 'laning region' which grows out from the wall 
indefinitely until the laning structure covers the entire 
range of the numerical grid. The value of W for laning 
states may thus be operationally defined as the density 
difference between the minina and maxima of the oscilla- 
tions at a distance sufficiently far removed from the wall. 
In practice, W was estimated from the variation of the 
profile around y = 150. The inset to FigUshows as a 
function of Pe for — 0.42, following the path indicated 
by an arrow in the main figure. For values of Pe close 
to, but above, the transition, W is well described by a 
square root. 

The phase diagram shown in FigU] is consistent with 
that calculated in Brownian dynamics simulations of 
charge stabilized colloidal dispersions (see Fig. 4 in j4l|). 
provided that the temperature used in the simulation 
study is identified as an inverse volume fraction. In 
[m temperatures both above and below the equilibrium 
freezing transition were considered and the nonequilib- 
rium order-disorder phase boundary found to vary con- 
tinuously through the freezing transition. In the present 
study we prefer to restrict ourselves to volume fractions 
below freezing (0^. = 0.494) in order to avoid the pos- 
sible complications which may arise from the presence 
of underlying mctastable states. A serious study of the 
complex interaction between crystal nucleation and ex- 
ternal flow goes beyond the scope of the present work. 

Finally, we note that the value Pccrit = 14.7 obtained 
from the present theory for = 0.45 is remarkably consis- 
tent with Brownian dynamics simulations performed at 
the same volume fraction (c.f. Figure. 3 in [36|). The sim- 
ulations predict that a layered structure emerges within 
the range Pe = 10 ^ 30. 
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FIG. 4: The phase boundary in the {<j) , Pe) plane separating 
the disordered phase from the ordered 'laning' phase (lines 
serve as a guide for the eye between calculated data points). 
The hard-sphere freezing transition at <^ = 0.494 is indicated 
by the broken line. The inset shows the order parameter W 
as a function of Pe for 4> = 0.42, following the path indicated 
by the blue arrow in the main panel. Above Pecrit the nu- 
merical data suggests a continuous transition with the order 
parameter varying aaW ^ (Pe— Pccrit)^ for small Pe— Pecrit- 



5. Bulk laning 

The results presented in the previous section indicate 
that the presence of the dynamical mean field term in 
(|17|) gives rise to an instability with respect to laning 
when the Peclet number exceeds a certain critical value. 
Although we have concentrated on the particular prob- 
lem of particles at a hard-wall, it would appear that 
the density inhomogeneities induced by the wall simply 
serve to 'seed' the generation of a laning structure for 
Pe > Pccrit- It may thus be anticipated that for su- 
percritical values of Pe, any kind of density fluctuation, 
regardless of its amplitude, will be sufficient to initiate 
laning. 

In order to test the above hypothesis we have solved 
P7p for a range of Pe numbers using the following initial 
guess for numerical iterative solution 

Pinit(y) = Pb + acxp{-b{y - yo'f), (33) 

for various values of the parameters a,b and yo. For 
Pe < Pccrit the perturbation is eroded during the it- 
eration procedure and yields the steady state solution 
pin) — Pb for all values of the parameter set (a, &, i/o)- 
For Pe > Pccrit any finite value of the parameter a 
is sufficient to seed the laning and a fully laned steady 
state solution, extending over the entire numerical sam- 
ple length, is obtained, regardless of the values of b and 
yo employed. In this sense, it would appear that, for 
supercritical states, any amount of 'numerical dirt' in 
the initial homogeneous density distribution is sufficient 
to generate a fully laned steady state. Moreover, we 
have confirmed that the values of Pecrit thus obtained 
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FIG. 5: Steady state sedimentation profiles of a relatively di- 
lute dispersion for Pe = 0, 31.416, 62.832, 103.673 and 125.664 
at a fixed value of the gravitational length ksT j glSra = 20. 
The black curve corresponds to Pe = and is obtained from a 
static DFT calculation at a fugacity z = 1.5. Due to the con- 
servation equation underlying our DDFT, steady state pro- 
files at finite Pe have the same adsorption as the profile for 
Pe = 0, i.e. particle number is conserved. Inset (a) shows 
the center-of-mass y (see Eg 1361 as a function of Pe. Inset (b) 
focuses on the region close to the wall and demonstrates fact 
that the contact value is independent of Pe. 



are entirely consistent with the phase boundary shown 
in FigHJ which was calculated in the presence of a hard- 
walL Given the above observations it would be of inter- 
est to perform a fully time-dependent solution of (|16p . 
Such a calculation, which we defer to a future publica- 
tion, would also enable predictions to be made regarding 
the timescale upon which lanes are formed and its de- 
pendence upon the supersaturation Pe — Pccmt- 



B. Influence of gravity 

We now consider adding an extra component to the 
external potential 



(34) 



where Am is the bouyant mass of a colloid and g is the 
gravitational field strength. We thus address the influ- 
ence of the shear flow ([TT|) upon colloidal sedimentation 
proflles. Choosing a fugacity z = 1.5 and gravitational 
length ^ = ksT/ gAm = 20 yields the equilibrium sed- 
imentation profllc shown in FiglSJ for which the local 
volume fraction remains rather low, even in the vicinity 
of the wall. As y increases, the local packing oscilla- 
tions give way to a monotonic decrease of the density. It 
may thus be anticipated that for flnite values of Pe, the 
flow kernel built into our mean-field theory will lead to a 
net transport of particles from regions of high density to 
regions of lower density until the scattering flux is bal- 
anced by the gravity-biased diffusion of particles towards 
smaller y values. 



The expectation of flow induced broadening of the sed- 
imentation proflles is conflrmed by the steady state re- 
sults shown in FigO Note that the particle number (i.e. 
area under each of the curves in FiglSJ is conserved and is 
independent of Pe. The canonical nature of the DDFT 
imposes that the broadening of the sedimentation pro- 
fllc with increasing Pe is accompanied by an overall de- 
crease in the local density within the range < y < 40. 
It is interesting to consider this change of the density 
distribution in view of the recently discussed violation 
of the fluctuation dissipation theorem (FDT) [4^ \^ in 
sheared systems. This violation was expressed in terms 
of the fluctuation dissipation ratio XpDT.f deflned as the 
ratio of response and thermal fluctuations for observable 
/ (see, e.g. [45| for details). In equilibrium, one has 
^FDT,f = 1, while under shear, XpDT.f < 1 is observed. 
Since, by deflnition, the ratio XpDT.f is proportional 
to (fc^T)"^ (when keeping response and fluctuations T- 
independent), one can also describe the FDT violation 
in terms of an effective temperature TpoT.t = T/Xpdt^{ 
which in turn is larger than T (Note that the de- 

pendence of ^FDT,f, and hence TFDT,f, on observable / 
is unclear.) Here, we are tempted to dcflne in analogy 
the centcr-of-mass ratio Xcom, 



Xcj 



(35) 



with 2/(7) the center-of-mass of the density distribution 
at shear rate 7, 



y 



J^dy yp{y) 
J^dypiy) ' 



(36) 



At low density, p{y) cx exp[— y/^] and y = ^ oc kgT, 
independent of shear. At higher density, p in Fig. [5] 
does not follow a simple exponential, but one still ex- 
pects y{Q) cx ksT, as long as packing effects are not too 
dominant. This conflrms the close analogy of our defln- 
tion of A'com to ATpDT.f : if for the system under shear, the 
ratio Xcom will be smaller than unity, one can formally 
interpret this in terms of an effective temperature larger 
than T. Inset (a) to Fig. 5 shows that the centcr-of-mass 
under shear is indeed larger then in equilibrium, i.e., we 
indeed have Xcom < 1 in accordance to the flndings for 
the ratio XpoT.f- The decrease of Xcom as function of 
shear rate resembles the behavior of XpD^ t, which was 
also found to decrease with shear rate [45,j47|. We fur- 
thermore expect that Xcom decreases with density and 



note Xn 



1 for — (where y{j) ^), as ob- 



served for XpY)T,f- The center-of-mass ratio hence shows 
the same overall properties as the fluctuation-dissipation- 
ratio. This suggests that both are driven by similar phys- 
ical processes. These flndings are also interesting in view 
of efforts towards a thermodynamic deflnition of an ef- 
fective temperature of the system under shear [4^. We 
realize that a comparison of the concrete values of Xcom 
and XpDT.f is not possible since the system under gravity 
is different from the bulk-systems studied in |45l447| . as 



10 




FIG. 6: Steady state sedimentation profiles of a dense dis- 
persion for Pe = 2, 4, 6 and 15 and at separations removed 
from tlie wall. The inset shows the local structure close to 
the wall. Although some local layering can be induced close 
to the wall at high shear rates, the gravitational force sup- 
presses the development of an extended layering phase. As in 
FiglH the contact values remains independent of Pe. 



the density depends both on y and 7. Future work on a 
single tagged heavy particle in a bath of density matched 
particles might prove more useful in this respect. 

Despite the broadening of the profiles as a function of 
Pe, the ultimate asymptotic decay can always be fitted 
by a Boltzmann decay, p{y — >■ c») ~ exp(— ?//^). This is 
expected since the density far away from the wall is low 
and gill) hence approaches the Boltzmann-distribution. 

The inset to Fig E] focuses on the region 0.5 < ?/ < 3. 
Despite the major changes in density distribution in- 
duced by the applied shear flow, it is striking that the 
contact value p[d/2) remains independent of Pe, in con- 
trast to our previous findings at g = (see Figs [2] and [3]). 
This is not surprising. For any finite g, the steady state 
density profile has an adsorption F = J^dy p{y) corre- 
sponding to the average number of particles in a column 
in y-direction with unit area in the x^-plane. In a gravi- 
tational field this column of particles thus exerts a force 
Tg/^m on the wall and determines the contact value of 
the density distribution. As particle number is conserved 
within the DDFT, it follows that the contact value of the 
density at the wall will be independent of Pe, as observed 
in our numerical solutions. The fact that the variation 
in contact value as a function of Pe is distictly different 
for the two cases Pe = and Pe 7^ is related to the 
singular behaviour of p{y — > 00) in the limit g ^ Q. 

Fig[6] shows further sedimentation profiles for a state 
with ^ = 10 and larger local volume fraction than those 
shown in FiglS] As previously, the profiles broaden with 
increasing Pe. Due to high local density in the vicinity 
of the wall, it may be expected that the layering transi- 
tion identified in our calculations at g = may become 
relevant at sufficiently high Pe values. The profile at the 
highest Pe value shown in FiglB] [Pe ~ 15) does indeed 



show the development of a layering structure close to the 
wall, similar to that in Figl^l However, the gravitational 
force acting on the particles supresses the development 
of long range oscillations and disrupts the formation of 
a macroscopic layering phase at any finite value of Pe. 
As for the profiles shown in FigEl the contact value at 
the wall remains independent of Pe over the range of 
parameters investigated. 



VIII. DISCUSSION 

We have applied dynamical density functional theory 
to calculate the density profiles of a colloidal liquid at a 
wall under shear flow. The chosen flow geometry served 
to highlight failings of the existing DDFT approach to 
driven states and a semi-empirical correction was pro- 
posed to reintroduce the missing physical mechanism. 
Calculations performed at various volume fractions and 
Peclet numbers have revealed that the new approxima- 
tion captures a non-equilibrium phase transition to an 
ordered laning state, for shear rates above a critical value 
of Pe. Moreover, sedimentation profiles are dramatically 
altered by the application of shear flow, which leads to an 
increase in height of the colloidal center-of-mass with in- 
creasing shear rate. The behavior of the center-of-mass 
ratio Xcora is in qualitative agreement with the previ- 
ously studied fiuctuation-dissipation-ratio under shear. 
The study of a single tagged heavy particle in a bath of 
density matched particles might be an interesting variant 
for the future. 

The mean-field correction to the advection term is 
presently rather empirical in character, arrived at using 
physical arguments, and it would be desirable to place 
this on a more rigorous basis, ideally as part of a system- 
atic scheme for improving the DDFT under external flow. 
Whether this is possible remains to be seen. In some 
sense, the present state of the theory for driven states 
is reminiscent of the early days of equilibrium DFT, for 
which the first attempts to go beyond the local density or 
square gradient approximation relied on the introduction 
of empirical weight functions to incorporate spatial non- 
locality 0. The insight gained from these studies proved 
very useful for the development of subsequent nonlocal 
approaches with a better foundation in statistical me- 
chanics. We thus hope that the present work may provide 
stimulus for further developments in applying DDFT to 
driven nonequilibrium states. 

The physical 'scattering' mechanism built into the 
present theory generates a nontrivial coupling between 
density inhomogeneities and external fiow, but has no 
effect on systems with a homogeneous density distribu- 
tion. While this is likely to be appropriate for certain 
colloidal systems, it may represent an approximation for 
others. Imposing shear flow on a homogeneous disper- 
sion generally leads to the development of flnitc nor- 
mal stresses, which have been associated with the phe- 
nomenon of shear-induced particle migration [4^ . It is 
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thus conceivable that a shear-induced drift of particles 
to regions of lower shear rate could result in a density 
gradient through the sample. Very recent experiments 
on PMMA hard-sphere-like colloids suggest that flow- 
concentration coupling can lead to a novel form of shear- 
banding [4^. However, the banding reported in [i^ only 
occurs for volume fractions above the glass transition, 
whereas the present work is focused purely on colloidal 
liquid states. 

Both the standard form of DDFT ^ and its advected 
extension ^ are based on an implicit adiabatic assump- 
tion which neglects the time taken for the (one-time) pair 
correlation functions to equilibrate following a change in 
the average density profile. It may thus be anticipated 
that in very dense systems, for which the structural a- 
rclaxation time becomes important, the pair correlation 
functions will be unable to keep up with changes in the 
density, leading to a breakdown of the adiabatic approx- 
imation. The fact that the structural relaxation time 
of driven dense states is determined by the inverse flow 
rate 7"^ (at least for states with < To) raises the 
interesting possibility that the adiabatic approximation 
may be more accurate when applied to calculate the re- 
sponse of dense systems to time-dependent changes in 
flow rate than to time-dependent changes in external po- 
tential. The present work has focused on steady state 
response and the next step in our research program will 
be to extend our studies to treat time-dependent shear 
flow. 

An important simplification of the present treatment 
is that hydrodynamic interactions have been neglected. 
This excludes from the outset the development of the 
hydrodynamically bound clusters which may form at 
very high shear rates and which have been s ugg ested 
as a possible mechanism for shear thickening [40[. As 
we focus here on the low and intermediate Peclet num- 
bers characteristic of the onset of shear thinning, this 
omission should not be too severe. More fundamental 
is the fact that the ordered phase observed in Brownian 



dynamics simulations |36j and captured by the present 
theory is apparently absent in Stokesian dynamics sim- 
ulations including the full solvent hydrodynamics [3]. 
It would thus appear that hydrodynamic interactions 
can render the ordered phase unstable. Nevertheless, 
we consider it important that any prospective theory 
of driven colloids be able to descibe first the simpler 
case of interacting Brownian particles, before seek- 
ing to refine this to include hydrodynamics at some 
level. While it may well be that the (approximate) 
incorporation of hydrodynamic interactions into the 
theory disrupts the laning behaviour reported here, we 
can at least be sure that such an improved theory has 
a sound physical basis and that the laning observed 
in Brownian dynamics fsof will indeed emerge should 
we choose to switch-off the hydrodynamics. Such a 
gradual theoretical development, adding new physi- 
cal aspects step-by-step, is important in developing 
a robust theory and tackling the fully hydrodynamic 
problem from the outset would be unlikely to deliver this. 
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